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Abstract —A mesh-dependent relation for the slip number in 
the Navier-slip with friction boundary condition for computations 
of impinging droplets with sharp interface methods is proposed. 
The relation is obtained as a function of Reynolds number, Weber 
number and the mesh size. The proposed relation is validated 
for several test cases by comparing the numerically obtained 
wetting diameter with the experimental results. Further, the 
computationally obtained maximum wetting diameter using the 
proposed slip relation is verified with the theoretical predictions. 
The relative error between the computationally obtained max¬ 
imum wetting diameter and the theoretical predictions is less 
than 10% for impinging droplet on a hydrophilic surface, and 
the error increases in the case of hydrophobic surface. 

Index Terms —Navier-slip, moving contact line, impinging 
droplet, finite elements, ALE approach 

1. Introduction 

Impinging droplets are encountered in many scientific and 
industrial applications such as spray cooling, inkjet printing, 
fuel injecting, etc. Simulating such fiows is complicated by the 
violation of the no-slip condition in the vicinity of the moving 
contact line, where the liquid-solid, solid-gas interfaces and 
the free surface intersect. The choice of the classical hydrody¬ 
namic “no-slip” boundary condition in the neighbourhood of 
the moving contact line leads to an unsatisfactory model that 
induce multivalued velocity field, refer ilD-ii- To alleviate 
this problem, often the contact line is allowed to move instead 
of imposing zero fiuid velocity at the contact line. A number 
of approaches have been proposed in the literature to move 
the contact line. In one of the approaches, the velocity of the 
moving contact line is prescribed as a function of the local 
dynamic contact angle 0, which is the angle between the 
liquid-solid interface and the free surface. Several models for 
the contact line velocity have been proposed in the literature, 
see Eggers et al. O for an overview. These models are mostly 
valid for wetting or perfectly wetting liquids. Eurther, the 
local dynamic contact angle is seldom available, and it varies 
for different fiow configurations. Therefore, this approach is 
hardly used in computations. Another approach is to allow the 
fiuid in the vicinity of the contact line to slip over the solid 
surface, refer ID, Q, 0 i.e., the relative velocity of the solid 
and liquid will be nonzero. To induce a slip, the slip with 
friction boundary condition 

(w - u) • Ty = • T(u,p) • Vy (I) 


is used, see for example, Gennes [O and Ganesan O • The slip 
boundary condition has first been proposed by Navier ifTOll . and 
later studied by Kundt et al. CD and Maxwell ca for gas 
dynamics. Here, (w — u) • Ts is relative velocity (tangential) 
of the solid and the liquid, and Ty • T(u,/7) • Vg is the shear 
stress of the liquid on the solid surface. Eurther, is the 
slip coefficient which defines the extent to which the no-slip 
boundary is relaxed. 

A relation between the Greenspan slip coefficient and the 
grid-spacing of the numerical scheme has been proposed by 
Moriarty et al. |[T^ for the moving contact line problem arising 
in dry wall coating. A number of theoretical and numerical 
investigations have been performed by several authors for the 
choice of the slip coefficient for specific moving contact line 
problems. Different expressions for the slip coefficient such 
as constants, functions of grid size, etc. have been proposed 
for specific moving contact line problems, refer |[T1, O, 
Q, ini-lED. Molecular dynamics simulations were often 
used to predict the slip coefficient for moving contact line 
problems, see m, El. In almost all of these simulations, 
the moving contact line is considered in channel fiows. Hence, 
the predicted slip values may not be generalized to all moving 
contact line problems, in particular, to impinging droplets. 
Even though the Navier-slip boundary condition ([T]) has been 
widely accepted for computations of moving contact line 
problems, a general mathematical expression or an empirical 
correlation for the slip coefficient is not available for imping¬ 
ing droplet simulations. The slip coefficient value need not be 
same for a droplet impinging on a same surface with different 
impact velocities. Often the slip coefficients for impinging 
droplets were identified on an ad hoc basis by comparing 
the numerical results with the experiments, see 0, El-El. 
The wetting diameter of the droplet has been used as a key 
parameter to identify the appropriate slip coefficient. A smaller 
value of the slip coefficient will reduce the wetting diameter, 
whereas a larger value increases the wetting diameter. Even 
though a deviation in the wetting diameter from the original 
value will induce a completely different fiow dynamics in the 
droplet, the equilibrium state of the droplet is not affected 
by the slip coefficient. However, an appropriate choice of the 
slip coefficient has to be used in computations in order to 
obtain physically accepted numerical predictions, especially, 
the dynamics of the fiuid fiow during the droplet deformation. 


It is the purpose of this paper to study the effect of the slip 
coefficient for different impact velocities and droplet sizes, and 
to compare the numerically obtained wetting diameter with 
experiments. Further, an expression for the slip coefficient 
is proposed. Apart from the choice of the slip coefficient, 
the inclusion of the contact angle into the model is very 
challenging. In particular, the choice of the contact angle 
value is very important in computations of impinging droplets, 
see Ganesan 1^ for a recent comparative study of different 
contact angle models. It has been observed that the equilibrium 
contact angle model is preferred for sharp interface methods. 
In discretization based numerical schemes (finite difference 
or finite volume or finite element methods), the contact angle 
is incorporated as a surface force, refer ||23l. Therefore, the 
measured dynamic contact angle need not be equal to the 
prescribed contact angle in the surface force until the droplet 
attains its equilibrium state. Consequently, the imbalance in 
the surface force induces a non-zero tangential velocity, and it 
necessitates slippage of liquid in the vicinity of the contact 
line. The above argument is another justification for the 
application of slip boundary condition in computations of 
moving contact line problems. 

An accurate approximation of the curvature and an appropri¬ 
ate discretization of pressure are essential to suppress spurious 
velocities in computations of free surface and interface flows, 
refer 1^ . In Eulerian approaches such as level set and 
volume-of fiuid methods, the free surface is not resolved by 
the computational mesh. Thus, an accurate calculation of the 
curvature and the conservation of mass are very challenging. 
Even though a separate surface mesh is used to explicitly 
represent the free surface in the front tracking method, the 
Navier-Stokes solver mesh does not resolve the free surface, 
and therefore the inclusion of the surface force is still chal¬ 
lenging. Alternatively, the free surface is resolved using the 
arbitrary Lagrangian-Eulerian (ALE) approach. Since the free 
surface is explicitly tracked in ALE approach, the surface force 
can accurately be incorporated in computations. Eurther, the 
inclusion of the contact angle is straight forward, refer ||23l. 
Even though handling the topological changes is very difficult 
in the ALE approach, it is possibly the most accurate approach 
for computations of free surface and two-phase flows when 
there is no topological change. Since the focus of this paper 
is to identify an appropriate expression for the slip coefficient, 
droplet impingement without any splashing and/or breakage 
is considered. Hence, the ALE approach is preferred in this 
study. 

The paper is organized as follows. The mathematical model 
and its dimensionless form of the governing equations are 
presented in Section 2. The used finite element scheme is 
briefly discussed in Section 3. The convergence study and an 
array of computations for impinging droplets are presented 
in Section 4. Eurther, a relation for the slip coefficient is 
derived and validated in this section. Einally, the findings are 
summarized in Section 5. 


II. Mathematical model 

We consider a spherical liquid droplet impinging on a 
horizontal surface, and the computation starts when the droplet 
comes into contact with the solid surface. Computations are 
performed until the prescribed time or until the droplet comes 
into the equilibrium after spreading and recoiling. A schematic 
representation of the computational model is presented in 
Eigure [B The liquid-solid interface and the free surface are 



Fig. 1. Computational model of a droplet impinging on a horizontal surface. 

represented by and Tp, respectively. Here, Oc denotes the 
contact angle, Tp, Vp are unit tangential and unit outward 
normal vectors on and T^, Vs are unit tangential and unit 
outward normal vectors on E5, respectively. 

A. Governing Equations 

The sequence of spreading and recoiling of an impinging 
liquid droplet is described by the time-dependent incom¬ 
pressible Navier-Stokes equations in a time-dependent domain 
Q.{t) C f G (0,7). 

^ + (u-V)u—— V• T(u,p) = f in fl(f) X (0,7) 

V-u = 0 in fl(f) X (0,7) 

( 2 ) 

where u denotes the velocity of the fiuid, p the pressure, p 
the density, 7 the given end time and f = (0,0,—g) the body 
force with gravitational constant g. The stress tensor T and the 
deformation tensor D for an incompressible Newtonian fiuid 
are given by 

T(u,p) := 2/iD(u) — pi, B(u) = ^ (Vu +Vu^), 
where p is the dynamic viscosity and I is the identity tensor. 

B. Initial and Boundary Conditions 

At time t = 0, we assume that the droplet is of spherical 
shape with diameter dp, and the initial velocity u(v,0) = (0,0,- 
Uimp{^))^ where Uimp is the impinging speed of the droplet. 
As mentioned in the introduction the Navier-slip with friction 
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boundary condition is imposed on the liquid-solid interface 
and it reads 

u-Vs = 0 on r 5 (f) X (0,/) 

Ts-T(u,/>)-vs = u-Ts on rs(r) X (0,/) 

The first condition is the no penetration boundary condition, 
i.e., the fiuid cannot penetrate an impermeable solid and thus 
the normal component of the velocity is zero. The second 
condition is the slip with friction boundary condition, i.e., on 
the liquid-solid interface, the tangential stress is proportional 
to the tangential velocity of the fiuid. Along the free surface, 
the force balancing condition 

T(u,/ 7 ) • Vp = Vr^ •Sri 7 on r/ 7 (r) x (0,7) 

is applied. Here, Vr^ and Vr^ • (•) denote the tangential gra¬ 
dient and tangential divergence, respectively, and are defined 
by 


and the Reynolds number, Froude number, Weber number and 
slip number are defined as 


Re = 


pUL 




1 

^ppU' 


III. Numerical scheme 


We use finite element method together with the ALE 
approach to solve the governing equations. We first derive 
a weak form of the Navier-Stokes equations. And then, we 
briefly describe the ALE formulation. After that, we discretize 
the weak problem in time and then in space. We briefiy present 
the numerical scheme here, and we refer to Ganesan et al. HU, 
ES-Ea for a detailed description. 


A. Weak formulation 

Let and be the usual Lebesgue and 

Sobolev spaces. We define the velocity space V and pressure 
space Q as follows : 


VrX-) = Pv^V(-), Vr^ ■(■) = & (PvA(-)), 

where = I — V/r (g) V 77 is the tangential projection. The 
surface stress tensor, Sri 7 1 ^ can be obtained as 

Sri7 = crPy^. 

Here, G is the surface tension. Further, the kinematic boundary 
condition 


u - Vf = w - V 77 on r77(r) X (0,7) 

holds, i.e. the normal component of the fiuid velocity on the 
free surface is equal to the normal component of the free 
surface velocity. 


C. Dimensionless form 

To write the Navier-Stokes equations in a dimensionless 
form, we introduce the scaling factors L and U as characteristic 
length and velocity, respectively. We define the dimensionless 
variables as 

V u ^ tU ~ lU p 

Using these dimensionless variables in the Navier-Stokes 
equations (O and boundary conditions and omitting the tilde 
after-wards, we obtain the equations in a dimensionless form 


|j^ + (u-V)u-V-T( u,j 3 ) = in Q(?) 

V - u = 0 in Q.{t) 

U-V5 = 0 onFsf) 

T5-T(u,;7)-V 5 = -PsU-Ts onr 5(0 

T{u,p)'Vf = — Vr^-Py^ onFfit) 

We 

u-Vi7 = w-Vp onFp{t) 


where the dimensionless stress tensor is given by 
T(u,p) = pD(u)-pI 


y = {v e //' (n(r))3 : v ■ = 0 on rs{t)} 

Q = {qeL^m))} 

To derive a weak form of the time-dependent incompressible 
Navier-Stokes equations, we multiply the momentum and 
mass balance equations by the test functions v G V and q G 
Q, respectively and integrate over After applying the 

Gaussian theorem for the stress tensor term and incorporating 
the boundary conditions, the weak form of the Navier-Stokes 
equations read: 


For given fl(0), u(x,0), find (u(x,t), p(x,t)) EVxQ such that 



+ a(u, u, v) -b{p,\)+b{q,u) = /(v) 


(3) 


for all V G V and q G Q. Here, 
2 

Ja{t) 


a{u,u,y) = — [ B(u):B(v)Jv+ [ (u-V)uvJv 

Re Ja{t) jQ{t) 

+ Pe {u-ts)(j-T:s)drs 

b(q,\) = / qV-ydx 

jQ-b) 

/W = f[ IPvF:VrfVfi;7 

FrJa(t) WeJrpit) 

+ — / cos{ec)\'Ts ds, 

We Jy, 


where jci denotes the contact line. We refer to Ganesan et 
al. 1241 for the inclusion of the contact angle. The contact 
angle model: 6c = 6e is used in all computations. The choice 
of equilibrium value in computations does not mean that the 
dynamic contact angle is fixed to the equilibrium value during 
the computations. Since the contact angle is included in the 
weak form as a natural boundary condition without imposing 
any condition on the geometry or on the contact-line velocity, 
the movement of the free surface in computations induces 
the hysteresis behaviour in the contact angle. A detailed 
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investigation on the effects of different contact angle models 
has been studied in Ganesan 1^ . and the equilibrium value 
is preferred for sharp interface methods. 


B. Arbitrary Lagrangian-Eulerian Approach 

Let At be a family of mappings, which at each t G [0,7) 
maps a point (ALE coordinate) Y of a reference domain Cl{t) 
onto the point (Eulerian coordinate) X of the current domain 
Q(r): 

A,iY)=X{Y,t) 


We assume that the mapping At is homeomorphic, i.e., At 
is invertible with continuous inverse. We also assume that 
the mapping is differentiable almost everywhere in [0,7). The 
reference domain Cl{t) can simply be the initial domain flo 
or the previous time-step domain when the deformation of 
the domain is large. Next, for a vector function u G C^(Q(t)) 
on the Eulerian frame, we define their corresponding function 
uecO(Q(r)) on the ALE frame as 

u:Q(t)—)^R, u:=uoAt, with u(Y,t) = u(At(V),t). 


Eurther, the time derivative of u on the ALE frame is defined 
as 


du 


! iQi(^) —^ M, 


5u 




r=A-^(x). 


We now apply the chain rule to the time derivative of uoA^ 
on the ALE frame to get 


(9u 

dt 


V 


(9u 

dt 


(X,t)- 


dx 

dt 


■ VvU = 


5u I 

dt \x 


W V;tU, 


where w is the domain velocity. Using the above relation, we 
write the Navier-Stokes equations in the ALE form as 


du I 


V • T{u,p) + ((u - w) • V)u = f, 


V u = 0. 


Since the free surface is resolved by the computational mesh 
in the ALE approach, the spurious velocities if any can be 
suppressed when the surface force is incorporated into the 
scheme accurately as discussed in Ganesan et al. Il26l . The 
application of ALE approach adds additional mesh velocity 
convective term in the model equations, and the mesh velocity 
needs to be computed at every time step. 


C. Axisymmetric formulation 

The computational domain of the considered problem is 
time-dependent and a very fine discretization, both in space 
and in time is needed to get an accurate solution. This 
requirement increases the computational costs in 3D. Since 
the considered domain is rotational symmetric, a 2D geometry 
with 3D-axisymmetric configuration is used. Thus, we rewrite 
the volume and surface integrals in @ into area and line 
integrals as described in Ganesan et al. ||28]| . It allows to 
use two-dimensional finite elements for velocity and pressure. 
Eurther, it reduces the computational complexity in mesh 
movement. 


D. Discretization in time and space 

Various time stepping methods have been proposed in the 
literature. The Euler schemes are of first order and the Crank- 
Nicolson is of second order but the latter is not strongly A- 
stable. Thus, we prefer the second order, strongly A-stable 
fractional-step scheme, refer ||29l, ||30l. Next to guarantee the 
stability and high accuracy we prefer the inf-sup stable finite 
elements of second order. We use triangular elements that ap¬ 
proximates the complex domains more accurately. One of the 
popular inf-sup stable finite elements used in computations is 
the Taylor-Hood element, i.e., continuous piecewise quadratic 
approximations for the velocity and continuous piecewise 
linear for pressure, and it is used in this paper. Eurther, a fixed 
point iteration is used to linearize the Navier-Stokes equations 
at every time step. Einally, the system of linear algebraic 
equations arising from the linearized Navier-Stokes equations 
is solved using UMEPACK (direct solver), refer ED. 

E. Mesh movement 

A linear elastic mesh update technique is used to handle the 
mesh movement. After solving the Navier-Stokes equations 
in each time step, the boundary displacement is calculated 
using the fluid velocity on the boundary. Using the boundary 
displacement as a Dirichlet boundary condition, the linear elas¬ 
tic equation is then solved for the inner points displacement. 
Einally, the mesh is moved with the computed displacement 
to get the next time step domain, see Ganesan et al. ll^ for 
more details. 


IV. Numerical Results 

In this section, we present the numerical results for an ax¬ 
isymmetric spherical liquid droplet impinging on a horizontal 
surface. We first perform a mesh convergence study in which 
we vary the number of points on the free surface. After that, we 
perform an array of simulations for glycerin and water droplets 
impinging on a glass surface with different impinging veloci¬ 
ties. The flow dynamics of the droplet depends on the surface 
characteristics, Reynolds, Weber, Eroude and the slip number. 
Among these numbers only the slip number is a numerical 
model parameter. Thus, the effect of the slip number on the 
flow dynamics of droplet for different impinging velocities 
and liquids are studied. The appropriate slip number for each 
test case is identified by comparing the numerically obtained 
dimensionless wetting diameter with their corresponding ex¬ 
perimental result presented in the literature. Based on the 
identified slip values, a correlation for the slip number in 
terms of the mesh size, the Reynolds and the Weber number 
is obtained. An array of simulations are performed by varying 
the equilibrium contact angle to check the applicability of 
the proposed slip relation for hydrophilic and hydrophobic 
surfaces. The maximum wetting diameter obtained from the 
simulations using the proposed slip relation are compared with 
the analytical values and other experiments to validate the 
proposed slip relation. 
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A. Mesh convergence study 

In this section we perform a mesh convergence study for 
the proposed numerical scheme. Space discretization is a very 
important aspect in CFD simulations in order to obtain accu¬ 
rate numerical results. Numerical simulation with extremely 
small mesh size is ideal to the continuum problem but it is not 
possible in practice due to the limited computational resources. 
We use open source package Triangle for mesh generation, 
which is based on constrained delaunay triangulation and the 
constraint we impose in our problem is the number of points 
used to track the free surface. In order to identify a feasible 
mesh size, we perform an array of simulations with a test 
example by varying the number of points on the free surface. 

We consider a spherical water droplet of diameter 
do = 2.7 mm. We take the characteristic length L = dof2 = ro, 
characteristic velocity U = Uimp and the dimensionless numbers 
used in the computations are Re = 1573, We = 25, Fr = 104 
and 6e = 75°. Five variants for the free surface points have 
been used which are as follows: (i) LO : 25, (ii) LI : 50, 
(hi) L2 : 100, (iv) L3 : 200 and (v) L4 : 400. First, we 
use a constant slip number (jSg = 30) in all the five variants. 
From Figure O we observe that the wetting diameter increases 
with increase in the number of points on the free surface. 
Hence, we cannot obtain convergence using a constant slip 
number. But from the wetting diameter curve, we can infer 
that the slip number has to be chosen in such a way that the 
wetting diameter is reduced with increase in the free surface 
points. Also, we know that the wetting diameter decreases with 
increase in the slip number value and the mesh size decreases 
with increase in the free surface points. Hence, we need to use 
a mesh-dependent slip number. Now, we perform computations 
using a mesh-dependent slip number, jSg = j8/ho, where ho is 
the initial size of the mesh onTp. For the values of slip number 
used in the computations, refer to Table [D The computationally 
obtained wetting diameter and the dynamic contact angle are 
shown in Figure [51 From Figure Oa), we can observe that 
there is almost no influence of the free surface points on the 
wetting diameter. As ho tends to zero, jSg tend to infinity which 
leads to the no-slip condition. Hence, the slip number can 
be interpreted as an artificial friction/slip introduced in place 
of no slip condition for moving contact line problems. From 
Figure Ob), we can observe that the free surface points have a 
significant infiuence on the dynamic contact angle. However, 
we can see convergence with L3 and L4 meshes. Since, our 
aim is to accurately capture the flow dynamics of the droplet, 
in all the subsequent computations we use L3 mesh, i.e. 200 
points on the free surface. 

B. Glycerin droplet 

In this section we consider glycerin droplets impinging 
perpendicularly on a smooth glass surface with equilibrium 
contact angle of 15°. The used values of physical param¬ 
eters are : p = 1220 kg m“^, p = 0.116 N s m~^ and 
(7 = 0.063 N m“^. Further, we take U = Uimp, L = dol2 = ro, 
Pe = jS/ho with ho = 0.01557859 and g = 9.8 m s“^. The 
impinging velocity of the droplet is varied between 1.41 m s“^ 


TABLE I 

Different cases oe eree sureace points used eor convergence 

STUDY ON A SPHERICAL LIQUID DROPLET 


Variant 

Points on F/r 

ho 

p 

R -1 

LO 

25 

0.12462872 

0.467343 

3.75 

LI 

50 

0.06231436 

0.467343 

7.5 

L2 

100 

0.03115718 

0.467343 

15 

L3 

200 

0.01557859 

0.467343 

30 

L4 

400 

0.007789295 

0.467343 

60 


TABLE II 

Dieeerent cases oe glycerin droplet used in this work 


Case 

Re 

We 

Fr 

Uimp(m s“^) 

Ps (identified) 

A 

18 

47 

166 

1.41 

2000 

B 

24 

81.5 

286 

1.854 

750 

C 

31.5 

140 

492 

2.43 

300 

D 

37.5 

201 

706 

2.91 

200 

E 

44.5 

285.5 

1002 

3A1 

125 

F 

61 

528 

1856 

4.72 

25 


and 4.72 m s“^. The obtained corresponding dimensionless 
numbers using the above parameters are given in Table [III 
Computations are performed till the dimensionless time 100 
with a time step length of 0.0005. For each case in Table [III nu¬ 
merical simulations are performed with different slip numbers. 
The formation of secondary droplets (topological changes) is 
not considered and it is the reason for the choice of this specific 
range of impinging velocity of glycerin droplets. 

We first study the infiuence of the slip number on the 
wetting diameter. Greater the value of slip number implies 
greater the effect of artificial friction. Hence, jSg -A oo im¬ 
plies no slip and jSg ^ 0 implies free slip. In Figure [H 
the dimensionless wetting diameter obtained with different 
slip numbers for the case F is in good agreement with the 
experimentally observed values till the dimensionless time 
t = 1, i.e., till the initial spreading phase of the droplet. During 
the initial spreading phase, the effect of slip number on the 
fiow dynamics is very minimal. However, after this initial 
phase different slip numbers induce different fiow dynamics. 
For droplets with low slip numbers, the frictional resistance is 
less and hence the spreading velocity is higher when compared 
to the droplets with high slip numbers. Higher the spreading 
velocity, greater is the kinetic energy of the droplet. Also the 
wetting diameter directly depends on the kinetic energy of 
the droplet. Therefore, the maximum wetting diameter will 
be greater for low slip numbers and it can clearly be seen in 
Figure [U 

The viscosity of glycerin is two orders higher than that of 
water. High viscosity of droplet induces a large resistant to 
the spreading and recoiling of droplet. Hence, the glycerin 
droplet deforms slowly on a smooth glass surface and it 
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Fig. 2. Computationally obtained dimensionless wetting diameter (a) and dynamic contact angle (b) with different points on the free surface using constant 
slip number (jSg = 30) for the cases in Table U 




Fig. 3. Computationally obtained dimensionless wetting diameter (a) and dynamic contact angle (b) with different points on the free surface using mesh 
dependent slip number (jSg) for the cases in Table U 



compared with experimental results. 


takes long time to attain its equilibrium wetting diameter. 
Generally, glycerin droplet does not rebound much due to 
high viscous dissipation. Also the equilibrium contact angle 
will influence whether the droplet will recoil or not after 
reaching the maximum wetting diameter. The recoiling effect 
is not observed in the all the considered cases because the 
equilibrium contact angle is very small, i.e., Oe = 15°. The 
maximum wetting diameter is same as the flnal equilibrium 
wetting diameter in whole range of the investigated impinging 
velocities. Also the maximum wetting diameter increases with 
increase in the impinging velocity of the droplet. We can 
observe in Figure lU that the slip numbers have a signiflcant 
influence on the flow dynamics of droplet after the initial 
spreading phase. Hence, choosing an appropriate value for 
slip number in the computations is very essential indeed. On 
comparing the numerical simulations with experimental results 
from Sikalo ll32l . we identifled an appropriate value for the slip 
number in each test case. The identifled values of slip number 


(Pe) are 2000, 750, 300, 200, 125 and 25 for the cases A, B, 
C, D, E and F, respectively, and are presented in Table [III Note 
that all the slip number (jSg) values indicated above are of the 
form Pe = jS/hp with ho = 0.01557859. We can also observe 
that the identifled values for the slip number decreases when 
the impact velocity increases for glycerin droplet. 

C. Water droplet 

In this section we consider a water droplet impinging 
perpendicularly on a smooth glass surface with equilibrium 
contact angle of 10°. The used values of physical parameters 
are: p = 996 kg m“^, p = 10“^ N s and a = 0.073 N m“^. 
The impinging velocity of the water droplet is varied between 
0.764 m s“^ and 2.96 m s“^ The corresponding dimensionless 
numbers obtained using the above parameters are given in 
Table jllll Computations are performed till the dimensionless 
time 10 with a time step length of 0.0005. For each case in 
Table [Till numerical simulations are performed with different 
slip numbers. Although the water droplet has comparable 
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TABLE III 

Different cases oe water droplet used in this work 


Case 

Re 

We 

Fr 


Ps (identified) 

G 

915 

9.5 

50 

0.764 

100 

H 

1573 

25 

104 

1.17 

30 

I 

1820 

38 

196 

1.52 

20 

J 

2810 

80.5 

330 

2.09 

10 

K 

2910 

97 

502 

2.429 

7 

L 

3545 

144 

746 

2.96 

4 


initial droplet diameter, equilibrium contact angle, surface 
tension and density to that of the glycerin droplet, its viscosity 
is two orders lower than that of glycerin. Due to its low 
viscosity, the droplet spreads more than that of glycerin. 
The rate at which water spreads is much higher compared 
to glycerin and this is the reason we have performed the 
computations only till dimensionless time t = 10. In certain 
cases, the computations are stopped due to the formation of 
secondary droplets (topological changes) or due to dry out of 
the droplet on Fs at the axis of symmetry. Because of low 
viscosity of water, we have chosen an even lesser range of 
impinging velocity for water droplet in this study in order to 
resist the early formation of secondary drops or the occurrence 
of splashing. 

The numerical result for the case J in Table |III| is shown 
in Figure HI During the initial spreading, we can observe a 
significant infiuence of the slip number on the fiow dynamics. 
This is in total contrast to what we observed in the glycerin 
droplet. This can be attributed to the fact that water spreads 
swiftly compared to glycerin because of significantly lower 
viscosity. For a given impinging velocity, the wetting diameter 
is higher for low slip numbers which was also the case with 
glycerin droplet. Also with increase in impinging velocity, 
the wetting diameter of the spreading droplet increases. The 
recoiling effect is not observed because of the choice of a 
small equilibrium contact angle, i.e., 6e = 10°. From Figured 
we observe that slip numbers have a significant infiuence 
on the fiow dynamics of the water droplet. On comparing 
the numerical simulations with the experimental results from 
Sikalo (3211 and Roux et al. (331, we identified an appropriate 
value for the slip number for each test case. The identified 
values of the slip number(j3e) are 100, 30, 20, 10, 7 and 4 for 
the cases G, H, I, J, K and L, respectively, and are presented in 
Table Hill We can also observe that the identified value for the 
slip number decreases when the impact velocity increases for 
water droplet which was also observed in glycerin droplet. On 
comparing the slip numbers for glycerin and water droplets 
with comparable impinging velocities, the slip numbers for 
glycerin droplets are almost two order higher than that of water 
droplet. Figure [5] depicts the magnitude of the velocity and the 
pressure contours of an impinging droplet (Case H in Table Hllb 
at dimensionless time instances t = 0.1, 1, 2, 5 and 10. 


D. Relation for the slip number 

Slip is a crucial factor in spreading of moving contact 
line problems. The numerical method introduces a slip at the 
discrete level, effectively introducing slip length on the order 
of the mesh size. Several authors na, d-ED have reported 
a convergence breakdown with the grid refinement and they 
overcame this by using a mesh-dependent slip for numerical 
solutions of moving contact line problems, which we observed 
in the earlier mesh convergence study. A relation between 
the Greenspan slip coefficient and the grid-spacing of the 
numerical scheme has been proposed by Moriarty et al. (T^ 
using curve fitting for the moving contact line problem arising 
in dry wall coating. Hence, this gives us the motivation to find 
a relation for the slip number applicable to impinging droplets. 
In the previous sections, we identified appropriate slip values 
for several test cases of glycerin and water droplet impinging 
on a glass surface. The dynamics of wetting for glycerin and 
water are not the same, e.g. different time scales for reaching 
maximum wetting diameter which is due to different viscosity 
in both liquids. However, the whole area of dynamic wetting 
has been motivated by developing models which are capable of 
describing widely varying wetting phenomena with the same 
set of parameters. Hence, this motivates us to obtain a relation 
for the slip number applicable to any liquid. 

We have studied the infiuence of the slip number on the fiow 
dynamics using the dimensionless wetting diameter which is 
also known as spread factor. The spreading behavior largely 
depends on the viscous and capillary forces of the droplet. 
The dimensionless numbers which account for these forces 
are the Reynolds and the Weber number, respectively. From 
the slip values, we observe that with increase in the Reynolds 
number, the slip number decreases and the decrease is quite 
rapid indicating that the relation may not be linear but could 
be exponential. The same behavior is also observed with the 
Weber number. Both the Reynolds and the Weber number 
play a major role in determining the spreading behavior. The 
dimensionless number which represents the relative effect of 
the viscous forces and surface tension is the capillary number, 
i.e., the ratio of Weber number to Reynolds number. However, 
trying to find a relation between capillary number and slip 
number will lead to the assumption that the relative effect of 
viscous forces and surface tension would be the same for all 
the droplets, which may not be true always. Hence, using the 
identified slip values for several test cases, we obtain a relation 
for slip number in terms of the mesh size, the Reynolds and the 
Weber number. For curve fitting, we used an online package 
called “Labfit”. Upon fitting, we have obtained the following 
relation. 

= J3 = aRe^ + AWe^ (4) 

ho 

where 

a = 4.796842276577 X 10^ 7 =-3.339370111853, 

X = 2.021796892969 x lO', 5 = -1.142224345078. 

Note that we have used L = ro in computations and the fit is 
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Fig. 5. Magnitude of the velocity (a) and the pressure (b) contours of a impinging droplet (Case H in Table HiH at dimensionless times t = 0.1, 1, 2, 5 and 
10 from the top. 
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TABLE IV 

Different cases oe equilibrium contact angles eor Water 
DROPLET WITH RE = 1573 


Oe 

WdN 

WdA 

Relative error (%) 

10 

3.8037 

4.0774 

6.71 

20 

3.6435 

4.0018 

8.95 

30 

3.4783 

3.8868 

10.51 

40 

3.3148 

3.7453 

11.49 

50 

3.1469 

3.5901 

12.34 

70 

2.7945 

3.2789 

14.77 

90 

2.5280 

3.0062 

15.91 

100 

2.4476 

2.8911 

15.34 

120 

2.3031 

2.7055 

14.87 

140 

2.1905 

2.5777 

15.02 


TABLE V 

Comparison oe numerical and experimental results eor 

VALIDATION 


Re 

We 

Fr 

Oe 

Pe 

WdE 

WdN 

Er 

1042 

29.5 

2257 

27 

27.19 

3.47 

3.45 

0.58 

1649 

59 

2846 

27 

12.32 

4.07 

4.07 

0 

2129 

85.5 

3163 

27 

8.06 

4.2 

4.39 

4.52 

2528.5 

109.5 

3342 

27 

6.08 

4.3 

4.6 

6.98 

1042 

29.5 

2257 

62 

27.19 

3.15 

2.91 

7.62 

1649 

59 

2846 

62 

12.32 

3.56 

3.54 

0.56 

2129 

85.5 

3163 

62 

8.06 

3.82 

3.89 

1.83 

2528.5 

109.5 

3342 

62 

6.08 

4.1 

4.1 

0 


using the Reynolds and Weber number which also are in terms 
of L = tq. However, in the literature authors have used L = Jo¬ 
in such cases, the slip number shall be used as : j8e=j3/2ho, 
where j3 is obtained from the proposed relation (|4]). 


E. Validation of the proposed slip relation 


In this section we perform an array of computations by 
varying the Reynolds number, Weber number and the equi¬ 
librium contact angle to validate the proposed relation for 
the slip number. To validate the relation for any hydrophilic 
surface, we study the influence of contact angle on the flow 
dynamics of impinging droplet. We consider the cases H and 
L in Table |III1 We perform computations for these two cases 
with the respective slip values as predicted by the proposed 
relation 0 and by varying the equilibrium contact using flve 
variants: (i) 10°, (ii) 20°, (hi) 30°, (iv) 40° and (v) 50°. From 
Figure [6l we can observe that the effect of contact angle on the 
flow dynamics is quite signiflcant for both the flows. However, 
we can predict the maximum dimensionless wetting diameter 
for flows with varying contact angles using the following 
analytical relation, refer 13411 . 


(We+12)WdA = 8 + w4 


We 

3(1 —cosO) +4 _ 

vRe 


(5) 


The maximum wetting diameter obtained numerically 


(WdN) from these simulations are compared with the values 
predicted by the analytical expression (WdA) in Table [IVl 
We have performed the simulations for wetting and partially 
wetting liquids. It has also been established that the mean 
error in predicting the maximum wetting diameter by the using 
the analytical expression is 5.09% with a standard deviation 
of 5.05%. For the case with equilibrium contact angle of 
10°, we have a relative error of 6.71%. However, this is the 
case we had obtained the slip number based on comparison 
with experiments. We assume that the experimental results are 
accurate and hence we have a error in the maximum wetting 
diameter predicted by analytical expression to be 6.71%. In 
this case, the analytical expression over-predicts when com¬ 
pared to experimental results. Even though the relative error 
in most of cases in Table EyI is more than 10%, due to over¬ 
prediction of the analytical expression we expect the relative 
error to be less than 10% for the cases with equilibrium contact 
angles 6e < 90°, as our calibration of slip number is based 
on the experiments. For hydrophobic and super-hydrophobic 
surfaces, i.e. for 0^ > 90°, the proposed relation may not be 
valid which could be a future scope for research. Hence, we 
can use the obtained correlation for the slip number values for 
droplet impinging on a hydrophilic surface. 

We have used experimental data from Sikalo |[32ll and Roux 
et al. 13^ to compare the numerical results and derive the 
relation for slip number. We now compare the numerical 
results obtained using the proposed slip relation 0 with some 
other experimental data provided in Ford et al. Ea The 
considered test cases are indicated in Table El Note that we 
have used ho = 0.01557859 in the computations and we have 
considered only droplet impinging on a hydrophilic surface. 
From the Table [Vl we can observe that the relative error (Er) 
in the maximum wetting diameter between the experimental 
and the numerical result is safely less than 10% for all cases. 
This further validates the proposed relation for the slip number 
for hydrophilic surfaces. 

V. Summary and Future Work 

In this paper we proposed a free surface mesh-dependent 
relation (|4]) for the slip number used in the Navier-slip with 
friction boundary condition on the liquid-solid interface for 
computations of liquid droplet impinging on a hydrophilic 
surface. An array of numerical simulations of liquid droplet 
impinging on a horizontal surface are presented in the pa¬ 
per. Finite element simulations are performed using arbitrary 
Lagrangian-Eulerian approach to study the effect of slip 
number on the flow dynamics of glycerin and water droplet 
impingement. Computations are performed for different im¬ 
pact velocities and droplet sizes. Appropriate value for the 
slip number in each test case is identifled by comparing 
the numerical results with experiments. Further, using the 
identifled slip numbers for the given Reynolds, Weber number 
and the mesh size, a relation is derived for the slip number. 
The proposed relation is then validated by comparing the 
computationally obtained maximum wetting diameter with the 
analytical predictions and other experiments. The proposed 
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Fig. 6. Computationally obtained dimensionless wetting diameter with different equilibrium contact angles for the cases H and L in Table m 


relation is more reliable for droplet impinging on a hydrophilic 
surface. Moreover, for droplet impinging on hydrophobic and 
super-hydrophobic surfaces, the same relation for slip number 
may not be appropriate. However, this could still give a good 
indication of the range of the slip number to be used in 
computations. Further research has to be done for the choice 
of exact slip number for droplet impinging on hydrophobic 
and super-hydrophobic surfaces. 
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